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Abstract 



We investigate the stability, the dynamical properties and melting of a 
two-dimensional (2D) Wigner crystal (WC) of classical Coulombic particles 
in a bi-layer structure. Compared to the single-layer WC, this system shows 
a rich phase diagram. Five different crystalline phases are stable; the energet- 
ically favoured structure can be tuned by changing either the inter-layer dis- 
tance or the particle density. Phase boundaries consist of both continuous and 
discontinuous transitions. We calculated the phonon excitations of the system 
within the harmonic approximation and we evaluated the melting tempera- 
ture of the bi-layer WC by use of a modified Lindemann criterion, appropriate 
to 2D systems. We minimized the harmonic free-energy of the system with 
respect to the lattice geometry at different values of temperature/inter-layer 
distance and we found no temperature-induced structural phase transition. 
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63.20.Dj - Phonon states and bands, normal modes, and phonon dispersion 
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INTRODUCTION 



Classical charged particles confined in a single, two-dimensional (2D) layer localize into 
a hexagonal lattice (Wigner crystal) for sufficiently large densities and low temperatures.0 
Such a single- layer Wigner crystal (SLWC) has been realized, e.g., on the surface of liquid 
helium.! Colloidal particles dissolved in water and placed between two glass plates is another 
example of an experimental system where classical particles exhibit Wigner crystallization.! 
Electrons in 2D semiconductor heterostructures behave like quantum particles, but a strong 
perpendicular magnetic field quenches the kinetic energy, leading the system toward the 
classical regime. The quest for the observation of such a Wigner crystal has been the object 
of very intense work over the last decades.@ The meltine transition of the classical SLWC, in 
particular, has attracted a large body of invest igationB since the proposal of a dislocation- 
mediated melting mechanism, leading to the prediction of a continuos melting transition.!!"! 

Recently, a new 2D system has attracted the attention of several groups, namely, the 
Wigner crystal in a bi-layer structure. One of the peculiarities of the bi-layer Wigner crystal 
(BLWC), compared to the SLWC, consists in the rich phase diagram; it has been pre- 
dicted that different crvstalline structures are stable in different ranges of inter-layer dis- 
tance/charge density.crEj 

In the present paper we address the phase diagram of such bi-layer structuresJil We 
consider a BLWC of Coulombic particles evenly distributed between the two layers. In a 
classical BLWC, the inter-particle interaction can be characterized by a unique dimensionless 



parameter rj = dyn/2, where d is the inter-layer distance and n is the total charge density. 
7] represents the ratio between the inter-particle interaction between the layer and within 
each layer. Thus, in the classical case, the Hamiltonian of the system is only a function of 
1], which therefore determines completely the phase diagram at T = 0. This is in contrast 
with the equivalent quantum problem, where d and n do not scale out. The search for the 
stable structure of a classical BLWC, at T = and as a function of rj, is made easier by the 
following considerations: 1) due to the long-range interaction, the two lattices which occupy 
the two layers (sub-lattices) are staggered to maximize the inter-particle distance. Each 
lattice site sits at the center of a cell in the opposite layer; 2) there are two trivial limiting 
cases: at r/ = the two sub-lattices reduce to a SLWC, which is known to crystallize in a 
hexagonal lattice (one-component hexagonal lattice). At the opposite limit of large t] the 
two sub-lattices are weakly coupled and, therefore, the stable structure is constituted by two 
staggered SLWC (staggered hexagonal lattice). By comparing the static energy of several 
lattices, we find that five different phases are energetically favoured in different ranges of 
7]. The five structures, in order of increasing 77, are a one-component hexagonal lattice (I), 
a staggered rectangular lattice (II), a staggered square lattice (III), a staggered rhombic 
lattice (IV), and a staggered hexagonal lattice (V). These phases evolve one into the other 
through both first and second order phase transitions. 

There exist already a number of investigations of the T = phase diagram of the 
classical BLWC in various systems.illlHliJli Some of the previous investigations of the present 
systemS'El did not identify all five phases. In Ref. 112, the bi-layer electron system which 
forms in a single wide quantum well above a critical densityB was studied. In this case, 
the transition from the single layer to the bi-layer (and, at higher densities, to a higher 
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number of layers) is of first order, with t] which jumps from to ~ 0.27; therefore, the low-?^ 
phases I and II were not investigated in Ref. 12. Theoretical investigations suggest that also 
bi-layer structures in the quantum regime possess a complex phase diagram. In Ref. |1^ a 
structural instability is found in the strong coupling regime, when the tunneling probability 
decreases below a critical threshold as a consequence of layer separation. The five phases 
described above have been predicted to exist in a bi-layer quantum Hall system.lll In this 
case, the phase diagram is even more complex, because a spin structure is associated with 
the lattice structure. Other bi-layer structures with an even more complex phase diagram 
can be imagined, for example, with different densities of the two layers.i'0 

Phonon excitations of some of the above systems have been partially investigated, within 
the harmonic approximation. Phonon frequencies have been used to evaluate the zero- 
point energy and the critical density for cold melting of the quantum BLWC0'E1 also in the 
presence of a magnetic field. 

In this paper we address the non-zero temperature phase diagram of a classical BLWC, 
which has not been investigated so far. With this aim, we have systematically investigated, 
within the harmonic approximation, the phonon excitations of each of the five structures. 
This allowed us not only to determine the range of structural stability and the behaviour of 
the acoustical and optical modes, which could be subject to experimental observations, but 
also allowed us to estimate the melting temperature of the crystal through the Lindemann 
criterion. The latter states that melting takes place when the mean square displacement of 
the crystallized charges exceeds a certain fraction of the lattice parameter. It should be noted 
that in 2D the mean square displacement diverges logarithmically with the crystal size. On 
the other hand, the relative mean square displacement between nearest neighbours (NN's) 
is a well defined quantity and its value at melting has been determined from simulations of 
2D crystals.^ Therefore, the latter quantity has been used in the present paper to estimate 
the melting temperature. Our results show that, by changing rj at constant temperature, 
one can pass through alternating regions of crystalline and liquid order. 

The overlapping range of stability of different lattice structures along the rj axis suggests 
the possibility that temperature-induced structural phase transitions take place, before the 
melting temperature is reached. To investigate this last possibility, we have minimized the 
harmonic free-energy of different structures at increasing temperatures at fixed 1], and we 
have found no evidence of such a temperature-induced structural phase transition. 

The paper is organized as follows. In Sec. | we investigate the zero temperature phase 
diagram. In Sec. |I| we calculate the phonon excitations of the systems. The non-zero 
temperature phase diagram is investigated in Sec. |T|. Results are discussed and summarized 
in the last section. 



I. ZERO TEMPERATURE PHASE DIAGRAM 

We consider a BLWC consisting of classical, spinless particles with total charge density 
n, evenly distributed over the two layers. The same form of the Coulombic interaction e^/r 
is assumed between particles in the same layer and in different layers. 

Electrons crystallized in the two layers constitute two sub-lattices which are equivalent 
by symmetry. When necessary, we denote the two sub-lattices by A and B. We consider 
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only the case of two layers of equal charge density = n/2 each; therefore, in the limit 
?7 — > 0, we recover a SLWC of density n. We consider the BLWC as a 2D lattice of N/2 
unitary cells, with two electrons per cell sitting on opposite layers. The primitive vectors are 
denoted by ai and a2, and the basis vectors are (0,0) and c; ai, a2, c, n^, and the vectors 
bi and b2 generating the reciprocal lattice, are listed in Table | for the five relevant phases. 
The equilibrium positions of the crystallized electrons in layer A and B are, respectively. 



Ra = iai + jaa 
Rb = isii + ja2 + c 



(1) 
(2) 



where i, j are integers. The total potential energy of the mobile charges due to the intra-layer 
and the inter-layer interaction is 



M+2 E 



1/2 



(3) 



The factor 1/2 accounts for double counting. Since the two layers are equivalent, and the 
origin can be chosen arbitrarily if we neglect surface effects, there are N/2 equivalent terms 
in each sum and the potential energy per particle E = Vp/N, therefore, reads 



E = -{Eo + Ej) 



where 



E() 



E 

R^O 



R 



(4) 



(5) 



represents the intra-layer interaction energy, and 



n [|R + c|^ + ci2' 



1/2 



(6) 



is the inter-layer interaction energy. Here R = isLi + ja2 and i? = |R| 
We follow Bonsall and Maradudin0 and rewrite as 



Eq = e lim 



E 

. R 



r-R| 



(7) 



where r is a 2D vector and r = |r|. It is convenient to define the following two functions 

giq-{r-R) 

R 



To (r, q) = e" 
Ti (r, q) = e 



|r — R| r ' 

giq-(r-R+c) 



R 



[|r-R + 



rf2 



1/2 • 



(8) 
(9) 
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from which Eq and Ej are obtained 



Eo = eMimTo(r,0), (10) 



r-»0 



Ei = e'Tj{0,0). (11) 

Due to the long range nature of the interaction, the lattice sums in Tq and Tj converge 
slowly. The Ewald technique is commonly used to overcome this difficulty, and consists in 
splitting the slowly convergent sum into two parts: the contribution of the first shells of 
neighbours is summed up in real space, while the contribution of the outer shells is summed 
up in reciprocal space. Both sums turn out to be rapidly convergent. The transformation 
to rapidly convergent sums over the real lattice vectors R and the reciprocal lattice vectors 
G is reported in Appendix 0. Here we state the final result 

To(r, q) = e-('^+«)-$ 0-^^^] + E ^''''''''^ |r - R 

+v^<l> {nrisr^) - i, (12) 

+y^Ee"''''^^"''^'^ (tt [ris |r - R + c|^ + 7]^]) . (13) 



R 

The functions $(x) and \l/(x, y), defined in Appendix 0, decay exponentially to zero for 
large x; therefore, Tq and T/ contain only rapidly convergent sums. 

The G = term in the first sum on the rhs of (|12]) and (|l^) give rise to a divergent term 
in Eq and in Ej, ensuing from the lack of charge neutrality. These terms are balanced by 
the interaction with a positive background, independently from the lattice geometry. Since 
the origin of the energy can be chosen arbitrarily, the divergent terms can be separated out 
and neglected in the calculation of the energy. This is done in Appendix ^ and the final 
result is 



Eo = 2e2y^ J2 ^ {-^^sR') - 2 , (14) 
Ej = e^v^ I E \^ + c\^ + V 

[ R 

+ E e'^"^ (t~' """^I + ^ erfcv^r] - C"' } I , (15) 



G^^O 



where G = |G| and erfc(a;) is the complementary error function defined in Appendix ^ 

We have calculated Eje^^fn (recall that n = 2ns) for the five different lattices listed 
in Table |. Phases I, III, and V are "rigid", meaning that, for a fixed density, the cell is 
uniquely determined. On the contrary, phases II and IV are "soft", because each of them 
contains a parameter, the ratio 02/01 of the length of the primitive vectors (aspect ratio) and 
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the angle between them, respectively, which can take on continuous values; at each value of 
T), this parameter is determined by energy minimization. 

Figure |1] shows the calculated energy per particle of phases I, III, and V as a function 
of rj; the energy of phases II and IV would be nearly indistinguishable on the scale of 
this figure and will be shown later in Fig. ^. As we discussed in the introduction, phase 
I is energetically favoured at very small rj, and it reduces to the SLWC at r/ = 0. At 
1] = we find E = —1.96052 e^\/n, which coincides with Ref. At the opposite limit 



of large r] the two sub-lattices become less and less coupled and the favoured geometry is 
composed of two staggered SLWC (phase V). Accordingly, the energy converges to the value 
Eq = —1.96052 e^-\/n/\/2, where the factor \/2 accounts for the reduced charge density. In 
the intermediate range of t] the energetically favoured structure is phase III. The geometry 
of the three phases is sketched in the insets of Fig. |I], close to the range where they have the 
lowest energy. 

Within the regions delimited by the open dots in Fig. |^, the energetically favoured 
structures are phases II and IV. These intermediate phases allow the lattice to pass from 
phase I to phase III, and from phase III to phase V, respectively. Note that phase II contains 
phase I and phase III as limiting cases, corresponding to the aspect ratios 02/01 = 
and a2/ai = 1, respectively. Analogously, phase IV contains phase III as limiting case 
for 6 = it/2. Therefore, the transitions I^II, II— *>III, and III— i>IV are continuous. Note, 
on the other hand, that there is no continuous way to pass from phase IV to phase V and, 
therefore, the transition IV— s>V is of the first order, i.e., dE/drj exhibits a jump. Apparently, 
this point has been overlooked in Ref. |12|, where the authors claim that the transition IV^V 
is continuous. The necessity for a first order phase transition when going from phase III to 
phase V has also been discussed, by use of general group theoretical arguments, by V. Fal'ko 
in Ref. |Tl]. 

Figure H(a) shows the transition I^II^III on an enlarged scale. Phase I is energetically 
favoured only in a very small range around rj = 0. As rj > 0.006, in fact, a rectangular unit 
cell with 02/01 < y/S (phase II) is energetically favoured. In the inset of Fig. ||(a) we show 
how the aspect ratio 02/01 evolves in a continuous way during the transition; phase I evolves 
into phase III through an anisotropic shrinking of the rectangular unit cell, and eventually 
02/01 = 1, corresponding to phase III, is reached at 77 = 0.262. 

The energy of phase IV is compared in Fig. ^(b) with the energy of phase III and V. For 
0.622 < Tj < 0.732 the staggered rhombic lattice (phase IV) has the lowest energy. As shown 
in the inset, increasing i] the angle 6 between the cell axes evolves continuously from 90°, 
corresponding to the square lattice (phase III), to 69.48°, and suddenly drops to 60°, which 
corresponds to phase V. The phase boundaries found above agree well with those found in 
Refs. and 0. 

To conclude this section, we give the asymptotic expressions of the static energy, for 
small and large t], 

E 1 



-1.96052 + 



2v/2 



iTiri - 0.600434 77^ + 2.86713 rt for small 77, (16) 



E 1.96052 V 

e ' tor large ?7, (17j 



e^V^ V2 2v^ 

where v = -3(73/2)1/^ and w = {87r^/V^y/^ in Eqs. (0) and reproduce the 
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correct energy within 2 % for 77 < 0.3 and within 0.2 % for 77 > 0.5, respectively. We stress 
that the above expressions are not fitting functions, but have been obtained by a series 
expansion of ( |I5| ) with respect to rj in the relevant range. In Eq. (0) the linear term, 
ensuing from the last term in Eq. (|15D, is the only odd order term in the Taylor expansion 
and all higher order terms are even; the coefficients involve sums over the direct and the 
reciprocal lattice, which have been calculated numerically for the lattice of phase I. The 
coefficients in Eq. (|T^) can be obtained analytically, once one realizes that, for large 1], only 
the first shell of G's needs to be retained in ([T5|); higher order terms are proportional to 
e~'"^ and decay faster for large 1]. In Ref. |T2| two fitting expressions for the classical energy 
were given; however, we found that none of them have the correct limiting behaviour. 



II. DYNAMICAL PROPERTIES 

In this section we calculate the frequencies of the phonon excitations of the five different 
phases within the harmonic approximation. 

For a general lattice, the square of the phonon frequencies are the eigenvalues of the 
dynamical matrix defined by0 



1/2 



(18) 



where R^^ is the position vector of the K-th particle in the l-th cell of the crystal, and m^; 
its mass. The quantities {Ik, I'k,') are the force constants defined by 



(f)a/3 (^/t, i'l^') = dad/3(f) (Rl^ - R^'k' 



(19) 



where (R/^ — Rz'k') is the two-body inter-particle potential. Here and in the following we 
use the notation 5a-F(x) = (9F(x')/9x'„|x=x', where Xa is the a-th component of the vector 
X. Due to translational invariance, the force constants satisfy the sum rule@ 



0,^(/ft:,/V) = 0. 



(20) 



Since each 2D unit cell of the BLWC contains two electrons, the dynamical matrix is a 
4x4 matrix which we write in block form as 



D 



(21) 



where D^^, D^-^, D-^^, and D-^-^ are 2x2 matrices. Applying (|I^) to the BLWC and using 
translational invariance, we obtain the matrix elements of D"^^ and 



j^AA 


;q); 


a/3 


1 


j^AB 


;q); 


a/3 


1 



^0„/3(R)e" 

R 



-iq-R 



(22) 
(23) 



R 
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where rup is the electron mass, and the force constants are 



and 



bc,p(R) = dad (3 — , 

ri 



()„/3(R- c) = dadp- 



R^O, 



iR-cr + d2 

Using Eq. (pOl), we find the force constant for R = 



1/2- 



(24) 



(25) 



6,/3(R = 0) = - 



^0«^(R) + ^0,/3(R-c) 

R^^O R 



(26) 



Furthermore, since the two sub-lattices are equivalent, we have D"^"^ = D^^ and, using 
Eq. (pD, D^^ = [D^^]t. 

It turns out to be convenient to define 



S^-^(q) 
S^^(q) 



-iq-R 



g-iq-(R-c) 



R 



iR-cr + ci2 



1/2 • 



which can be obtained from To and Tj 



S^^(q) 



= -e^ limdo^d/sTo (r, q) 

ap r— ►O 



Then the matrix elements of the dynamical matrix can be written 

S^^ (0) + S^^ (0) - S^^(q) 
-S^^(q)l . 





;q) 


1 








j^AB 


;q) 


1 






rrie 



(27) 
(28) 

(29) 
(30) 

(31) 
(32) 



Using the rapidly convergent form for Tq and Tj, as given in (12) and (|T3D, allows one to 
write down the matrix elements of S^"^ and S^"^ explicitly 

S^-^(q) 



a/3 



-E(q + G)jq + G),JM^ 



Txn^ 



(33) 



Rt^o 



S^^(q) 



E (q + G), (q + G), VP f e-«- 



R 



(34) 
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where we have defined 



V^p{X^) = d^dp^iX^) = 2™, {5,^$' (X^) + 27:nsX^Xp^" (x^)] . (35) 

In general D is a complex hermitian matrix. However, since in the BLWC the two sites of 
each cell are occupied by identical particles, it is possible to apply a unitary transformation 
which results in a real symmetric matrix.E3 If we denote with I2 the 2x2 identity matrix, 
the transformation 

results in 

D = UDU-^ = ( ^^i,^}^ ^ V^^T. ) , (37) 
\^ Re Dab D^a - Im D^b / ' ^ ^ 

where Re ^ab and Im T)ab are the real and imaginary parts of D^^. Note that Im D^^ = 
for a lattice with inversion symmetry. This applies to all phases, except for phase V. 
Finally, we solved the set of four linear equations 



(D(q) 



2 



l4)e(q,j)=0, (38) 



where I4 is the 4x4 identity matrix, cjq ,,■ is the frequency of the j'-th phonon mode (j = 
1, . . . ,4) with wavevector q, and e(q, j) its eigenvector. Equation (pHf) is equivalent to the 



diagonalization of the 4x4 matrix D, which provides the four eigenvalues cUq^- at each point 



q in reciprocal space. For a lattice to be stable, it is necessary that cUq^ > 0. 



Figure |^ shows the frequencies tUqj (or, when tUqj is negative, its imaginary part) for 
phases: (a) I, (b) III, and (c) V, and their evolution with 77. Frequencies are given in terms 
of the characteristic frequency uji = e^n'^'^^/me, which depends on the density and not on 
the lattice geometry. Phonon dispersions are shown along the high symmetry directions in 
reciprocal space. The high symmetry points are labeled according to the insets. We recall 
that in a SLWC the transverse acoustical (TA) and the longitudinal acoustical (LA) modes 
vanish at the F point as q and g^/^, respectively.^ Thus the sound velocity of the LA mode 
is infinite. The latter behaviour is a general property of a 2D Coulomb plasma^ and does 
not depend on the lattice geometry, nor on d, as is clear from a comparison of the three 
panels in Fig. |^. The remaining two (optical) modes ensuing from (|38| ) are peculiar to the 
BLWC and correspond to out-of-phase vibrations of electrons in opposite layers. 

Starting from the top panel, it is shown in Fig. |^ that, as 77 is increased, the TA mode of 
phase I softens until, above a critical value of 77, the frequency becomes imaginary, indicating 
a lattice instability. For r] between 0.262 and 0.622 phase III (square lattice) is energetically 
favoured, according to Fig. |^. We recall that in a SLWC the square lattice cannot exists, 
since it has an imaginary TA branch.^! Figure ^(b) shows, on the other hand, that in the 
BLWC the square lattice is stable for a certain range of rj. For t] below a critical value, 
however, the TA branch softens along the FM direction and eventually becomes imaginary, 
as is expected from the fact that in this limit the BLWC tends to the SLWC; at the opposite 
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limit of large t], the TA branch softens along the FX direction and eventually phase III 
becomes unstable, as the BLWC tends to two separated SLWC's. 

Phase V is stable for large rj [Fig. 0(c)] . In the limit ^ oo we have the phonon 
dispersion curves of two uncoupled SLWCs; therefore, in Fig. |], each curve in the i] = oo 
case is doubly degenerate and all modes approach zero for q 0. For smaller rj, optical 
modes with a finite frequency at the F point appear; at the same time, the TA mode becomes 
softer and, eventually, becomes imaginary at ~ 0.6. 

The sound velocity of the TA mode, vta = d^^TA/ dq\q=o, along the in-plane directions 
(1, 0) and (1, 1), and for the five phases, is shown in detail in Fig. ^. The labels on top of the 
figure indicate the energetically favoured phase in each range of r], according to Figs. and |^; 
phase I is favoured only in a very small range around t] = and, therefore, is not indicated. 
The vanishing of the low-frequency modes in certain directions, shown in Fig. ^, sets a limit 
to the range of stability of each phase. Note that phases II and III have a soft mode at 
a value of rj which coincides with the value where the transition between the two phases 
takes place (the vertical dotted line). The same happens for phases III and IV. Therefore, 
the range where phase III is energetically favoured coincides with the range of stability of 
this phase. This has profound implications in determining the DLWC phase diagram at 
T 7^ 0, as we will show in the next section. Note also that in the range of rj where phase II 
is energetically favoured, both phases I and II are stable, i.e., they do not have imaginary 
phonon frequencies. Analogously, in the range where phase IV is energetically favoured, 
both phases IV and V are stable. The T=0 phase diagram, deduced from Figs. |I| and |^, 
and the range of stability, deduced from the softening of phonon modes, are summarized in 
Fig. |. 

Figure ^ shows the evolution of the optical frequencies at the F point, Uopt, with 77. 
For phases III and V the two optical frequencies are degenerate. It has been noted that the 
detection of the exponential decay of the optical modes at large rj could serve as a fingerprint 
of the solid phase in a bi-layer structure.liil Note also the different behaviour of the optical 
modes between phase I and II, and between IV and V, which may be used experimentally 
to distinguish between the different possible phases. 

The dependence of the sound velocity and of the optical modes at the F point upon rj 
has been fitted to simple analytical expressions in the low t] (below ~ 0.2) and in the large 
T] (above ~ 0.7) range. The sound velocity vta (see Fig. 0) has been fitted to 

vtaV^/ uJi=po+ P2V'^ + PiV^ (39) 

for small rj, and to 

VTAy/n/^^i =Po+ P2V~'^ + PaV'"^ (40) 

for large rj. The coefficients pi for phases I, II (small rj) and V (large rj) are reported in 
Table |n|. The frequencies Uopt (see Fig. H) have been fitted to 

^opt/^i = go + q2ri^ + Qif]'^ (41) 

for small 77, and to 

uJoptl^i = qie-''" (42) 
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for large rj. For 7] ~ 0.262, close to the boundary between phase II and phase III, the optical 
modes of phase II have a singular behaviour. In this range we have fitted Uopt to 



i^opt/uJi = qo + gi(g2 - vy^- (43) 

The coefficients are reported in Table |T|. The agreement with the full calculations did 
not improve by adding odd powers of rj in Eqs. (P5|), (plD, and (^Tj). All the above fitting 
functions give the correct values with an accuracy better than 0.6 % in the relevant ranges 
of 1]. 

In Fig. 1^ we report the evolution with rj of the phonon density of state (DOS). At each 
value of rj, we show the DOS of the phase which is energetically favoured at that value. 
At ?7 = and r] = oo the energetically favoured lattice are, respectively, the SLWC lattice 
with density n (phase I), and two uncoupled SLWCs with density n/2 (phase V); therefore, 
the corresponding DOS curves are equal up to a factor 2^/^ in the frequency scale. Note 
in Fig. 1^ the peak of optical frequencies which narrows at rj ~ 0.5, corresponding to the 
range of t] where the in-plane component of the average interaction of one particle with its 
NN's in the same layer and in the opposite layer are similar. Also note the low-frequency 
peak which moves to very low frequencies around t] ~ 0.3 and t] ~ 0.7. This behaviour is 
reminiscent of the softening of the TA mode of the square lattice (phase III) discussed above. 
The resulting high density of low-frequency modes suggests that very large fluctuations of 
particles around their equilibrium lattice sites are possible; correspondingly, a low melting 



temperature is expected in proximity of these points, as will be discussed in Sec. Ill 



III. PHASE DIAGRAM AND MELTING 

In this section we will be concerned with the non-zero temperature properties of the 
BLWC First, we will use the calculated phonon excitation frequencies to estimate the melt- 
ing temperature Tm via the Lindemann criterion. In principle, only order of magnitude 
estimates of Tm are expected from an harmonic theory, since anharmonic terms of the po- 
tential become important when crystal vibrations are so large that the lattice is near to 
dissolve. In the case of the SLWC, apart from simulations, analytical methods have been 
successfully used to calculate Tm by including anharmonic effects.E^J^ These methods assume 
that melting proceeds through the dislocation mediated mechanism proposed by Kosterlitz 
and Thouless,0 Halperin and Nelson,i and Youn^ (KTHNY theory). The ingredient of these 
calculations is the sound velocity of the TA mode of the lattice, which is assumed isotropic 
in the KTHNY theory and which is indeed the case in the simple hexagonal lattice of the 
SLWC, but not in the BLWC, where the TA mode, in general, is anisotropic, as is clear 
from Fig. ^. Therefore, in this work we will rely on the simple Lindemann criterion. We 
shall see that, taking the Lindemann parameter 6, defined in Eq. ( |4^ ) below, from existing 
simulations, effectively includes anharmonic effects into the theory to some extent. 

The Lindemann criterion states that, in a lattice of density n, melting occurs when 

^ = *^ (44) 
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i.e., when the mean square displacement of a lattice site around its equilibrium position (m^) 
exceeds a certain fraction of the mean inter-particle distance tq = l/y^vm. The brackets 
(. . .) represent the thermodynamic average; in our case the latter will be calculated within 
the harmonic theory. The parameter 6 is an input to the criterion, to be obtained from 



simulations or from some analytic theory. Equation P^ ) has been verified in simulations of 
several 3D systems.!! It is known, however, that (u^) is logarithmically divergent in 2D. On 
the other hand, the relative mean square displacement (|u(R) — u(R-|- a)|^), where u(R) 
and u(R-|-a) are the displacement vectors at lattice site R and at the NN site R + a, where 
a is the vector joining two NN's, is finite. Correspondingly, a modified Lindemann criterion 
can be defined 



(|u(R) -u(R + a)|') 



(45) 



The value of at melting has been calculated in simulations of melting in a SLWC and 
turned out to be ~ O.l.li^In principle, 5m may depend on the lattice geometry and the 
nature of the interaction; however, the Lindemann parameter has been found to be quite 
independent from the form of the interaction both in 2D0 and in 3D@ systems; therefore, 
we take 6"^ = 0.1, and independent from t] and from the lattice geometry. Small variations 
of 6m would not change our results qualitatively. 

The correlation function (|u(R) — u(R-|-a)|^) is calculated within the harmonic the- 
ory.0 Each lattice site in the BLWC has two types of NN's, in general at a different distance, 
and the number and distance of the NN's changes in a continuous way with rj. Accordingly, 
we calculate separately two (in general different) correlation functions 



Ml 



oi=x,y m=l...A/i 



AksT „ [e^(q,j)]' + [e;^(q,j)]' 



NnieMi 



E 



00, 



E 



sm 



q Rn 



(46) 



m=l...Mi 



M. 



E E (^f(o)-<(^) ) 



2 a=x,y m=l...M2 

ksT ^ 1 



^ {l-2[ef(q,j>f(q,j) + e^(q,J>f(q,j)]cosq-R„}, (47) 



m=l...A/i 



where fcs is the Boltzmann constant, u^*^^^ is the a-th component of the displacement vector 
in layer A(B) calculated at the origin (0) or at the position of the m-th NN. e^'''^''(q, j) is the 
a-th component of the eigenvector of the j-th mode, at point q, relative to the sub-lattice 
in layer A(B), R^ is the relative lattice vector connecting one site to its m-th NN in the 
same (Li) or in the opposite (L2) layer, and the sums over m are extended to the Mi (M2) 
NN's in the same (opposite) sub-lattice. 

Now we consider two limiting cases. For = 0, (|u(R) — u(R -|- a)|^) = L1 + L2, since all 
NN's are equivalent. At the opposite limit, 77 — 00, 

(|u(R) -u(R + a)n = 

Li, since the 

dynamics in one layer is not influenced by the sub-lattice on the opposite layers. Therefore, 
we write in general 
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(|u(R) - u(R + a)|') = Li + f{r])L2, (48) 
where the function /(r^) satisfies 

/(0) = 1, /(oo) = 0. (49) 

As f{r]) represents the infiuence of the oscillation in one layer on the oscillations in the 
opposite layer, we take f{ri) proportional to the in-plane component of the Coulombic force 
between two NN sites sitting in opposite layers. This is 

where c = |c|. Taking /(r/) proportional to -F|j((i), and imposing the conditions (^), we have 

m = ^ ,,3/, , (51) 

(1 + apT]^) ' 

where ap = (nc^) ^ is a dimensionless geometric factor which can be calculated from Table |. 

Inserting Eq. ( ^S]) in ( ^5] ) and using Tq = l/irng, we have calculated the melting temper- 
ature Tm, which is reported in Fig. || for the five phases. For the "soft" phases II and IV, 
Tm was calculated taking the T = value of the aspect ratio and 9, respectively. This will 
be justified later. 

In the studies of melting of the SLWC, the melting temperature is usually given in terms 
of the dimensionless parameter F^^ = e^y/wn/ KTm (the inverse of the vertical units in 
Fig. P), the ratio between the average Coulombic potential energy and the average kinetic 
energy. Experimentsi give Fj\,/ ~ 131 and simulationsil give Tm — 128. Using the harmonic 
value of the sound velocity at T = 0, the KTHNY theory gives F^ — 79. Our calculation, 
which is performed within the harmonic approximation, but uses 5m taken from simulations 
which, of course, include anharmonic effects, gives kBTM/e^\/T^ = 0.00925 at r/ = 0, cor- 
responding to Tm = 108. Therefore, our calculation, although overestimates Tm, partially 
includes anharmonic effects. In a full anharmonic theory, Li and L2, which in the harmonic 
approximation scale linearly with T, would increase more rapidly, especially close to the 
melting transition. 

Fig. § shows that the melting temperature has an oscillating behaviour as a function of 
rj. This is a consequence of the vanishing of the TA phonon modes at the phase boundaries 
II/III and III/IV, as discussed in Sec. |I|. Therefore, for fixed T 7^ and as function of rj we 
observe that alternating solid and liquid phases are possible, and the reentrant solid phase 
has a different lattice geometry each time. Furthermore, note from the inset of Fig. ^that, 
for large values of f], Tm approaches the melting temperature of a SLWC of density n/2 
from below. In certain experimental realizations of the BLWC it could be easier to change 77 
through a change in the charge density, keeping d constant. Therefore, in Fig. ^we show the 
calculated melting temperature in units of ksTMd/e^y/Ti. Note that in the classical regime 
the phase diagram is determined by t] and a dimensionless temperature, either kBTM/e^y/T^ 
or ksTMd/e'^y/Ti. In the quantum regime, instead, the kinetic energy term depends on 
the density alone and, therefore, the phase diagram must be drawn explicitly in the tree- 
parameter space {d,n,T). 
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The presence of different lattice geometries which are stable within the same range of rj 
suggests the possibility that, increasing T at fixed 77, the BLWC undergoes a structural phase 
transition, and, eventually, melts at a temperature appropriate to the high temperature 
phase. For example, it seems possible that for rj < 0.262 the BLWC evolves from phase 
II (with some value of the aspect ratio which minimizes the static energy at T = 0) to 
phase I (aspect ratio v^), as T exceeds some critical value, and eventually melts at a Tm 
appropriate for phase I. To investigate such possibility we have minimized the free energy 
with respect to the lattice geometry at fixed rj and T. The harmonic approximation of the 
free-energy in the high-temperature limit is 

i^(0 = ^(0 + fcBTElog%^, (52) 

where ^ is a parameter which defines a distortion of the lattice. There are two ranges of r) 
where more than one phase is stable with respect to lattice vibrations (see Figs. ^ and 
in the range 0.006 < r] < 0.262 phase II is energetically favoured, but also phase I is stable 
throughout this range. Therefore, in this range, we minimize F with respect to ^ = 02/01. 
In the range 0.622 < t] < 0.732 phase IV is energetically favoured, but also phase V is 
stable; therefore, in this range we take ^ = 6. Integration over reciprocal space in (|52D 
was performed numerically. We found that in both ranges the value of ^ which minimizes 
F is practically independent of the temperature and, therefore, coincides with the T = 
value. In other words, the phase boundaries between the different geometries in Fig. ^ are 
represented by vertical lines. Moreover, this justifies the fact that, in order to calculate Tm 
for the "soft" phases II and IV, we have used the T = value for the aspect ratio and 6, 
respectively. 



CONCLUSIONS 

The phase diagram of a classical BLWC, both at T = and at T 7^ was investigated, 
within an harmonic approach, by use of the Lindemann criterion and minimization of the 
harmonic free-energy. Five different crystalline geometries are stable in different ranges of 
inter-layer distance /charge densities. Moreover, at T = the five phases evolve one into 
the other through both continuous and discontinuous transitions. At T 7^ 0, alternating 
solid and liquid phases are possible, as one sweeps the inter-layer distance or the charge 
density. In particular, regions of hquid phase separate phase II from phase III, and phase 
III from phase IV. This has been shown to be a consequence of lattice instabilities induced 
by the vanishing of phonon modes at the phase boundaries. On the other hand, a first order 
transition line separates IV from phase V. 

An additional intricacy of the phase diagram in the small rj range has been pointed out 
by Vil'k and Monarkha0. In this limit the hamiltonian (|^) was mapped into the hamilto- 
nian of a binary mixture of particles sitting on a triangular lattice and interacting through 
a dipole potential. Therefore there is a possibility that a disordered phase appears, as the 
temperature is increased. They find two phase transitions. The low tempeature (ordered) 
phase is equivalent to our phase I. Above a critical temperature Ti the lattice can be seen 
as composed of three inter-penetrating triangular lattices, two of which are ordered and one 
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is disordered. Above a second critical temperature T2 the lattice becomes completely disor- 
dered. Of course the order-disorder transition vanishes as rj ^ 0, where the two sublattices 
become equivalent. 

A bi-layer electron gas can easily be realized in semiconductor heterostructures. Ili'§ Al- 
though our results have been obtained for a classical system, they can give some indications 
on the phase diagram in quantum bi-layer structures, provided that temperature fluctua- 
tions are interpreted as quantum fluctuations. Very recently, in Ref. ^ a re-entrant phase 
around rj ~ 2.6, analogous to ours in Fig. ||, was predicted in the {ri,rs) phase diagram, 
where is the dimensionless inverse electron density. Furthermore, our analysis of the 
phonon excitations and the analytical fitting that we have developed retain their validity in 
the quantum regime. 

In principle, the harmonic approximation used throughout this work is expected to fail 
when the temperature approaches the melting temperature. However, we have shown that, 
in the 77 = case, we obtain Tm in reasonable agreement with numerical simulations and ex- 
periments on the SLWC; therefore, we believe that inclusion of the anharmonic effects would 
not change our results qualitatively, as far as the melting temperature is concerned. We be- 
lieve also that the approximation of a structure independent parameter 6m in the Lindemann 
criterion should not change the nature of our findings. The harmonic approach could be a 
more severe approximation in the calculation of the free-energy and our investigation, there- 
fore, does not rule out completely the possibihty of temperature-induced structural phase 
transitions. 
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APPENDIX A: RAPIDLY CONVERGENT FORM OF Tq AND Tj 

The slowly convergent sums over lattice sites appearing in the definition of Tq and T/ 
[Eqs. (P) and (^)] cannot be used in a numerical calculation. Therefore, they will be converted 
into a rapidly convergent form using a generalization of the Ewald method.0 Formally, we 
proceed as follows. First, each term in the sum is decomposed in two terms, using the 
identity 
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r 



where 




(A2) 
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is the error function, erfc(x) = 1 — erf(x) is the complementary error function, and e is an 
arbitrary constant. The reason why we do so, is that erfc(x) vanishes exponentially for large 
values of the argument and, consequently, the lattice sum with this function as argument is 
sufficiently rapidly convergent. Then, the other lattice sum with argument erf(x) is mapped 
onto a sum over the reciprocal lattice, using the 2D theta-function transformation.0 
Using ([Al|) and the definition of Tq [Eq. (H)], we obtain 

jq-(r-R) 

% = e-'^ '^ E ——— erf {e |r - R|) 

g-'jq-R erfcfer) 1 , , ^, 

+ E I ^erfc(£|r-R|) + ^ - -. (A3) 

R^o |r - R| r r 

To convert the first sum on the rhs of ( |A3| ) into a rapidly convergent form, we substitute 
C, = t/X in the definition of the error function, which results into 

^^'(^^^-A Te-A^.^.^ (A4) 



A ^/7^ Jo 



with A = |r — R|. We plug ([A4|) into the first sum of ( |A3|) and we bring the sum under the 
integral. Next, we apply the 2D theta-function transformation^ 



J2 g-lr-Rl^g-iq-R. = !^ J2 e-li+GJ|'/4«'e"*('i+^)-^ (A5) 

R ^ G 

and the substitution t = |q + G| /2^, which transforms the first term on the rhs of ([A3[) into 



2nn. E e-('>^^)-^^^^^±^iiM. (A6) 



G 



|q+G| 



The final step is to choose a reasonable value for e, such that the lattice sums have a sufficient 
rapid numerical convergency. A convenient choice is e = rg ^ = i/vmj. Defining the function 

^{^) = y^erfc (v^) (AT) 

to simplify the final expression, we finally obtain Eq. ([T^). 

For Tj we proceed in a similar way. Let A^ = |r-R + cp + d'^; using (K]]), and the 
definition of Tj [Eq. (||)], we have 

Tj = e-''^ " ^'"^ [ erf (^A) + erfc {eX)] . (A8) 

R ^ 

Using the identity (^), the 2D theta-function transformation (^), and the substitution 
t = |q + G| /2^, the first term in ( [A8| ) becomes 

4 g-iG-c -i(q+G)-r „oo „, ^,2, 9 9 

V^r ^ |q+G| J|q+G|/2e 
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The integral can be performed analytically; using 



a 



e " erfc { x -\ — I + 

X 



erfc ( X — 



inserting e = and defining the function 



1 fW y 



m{x,y) = -^- [ev^erfc(v^ + yy) +e-v^erfc(v^- 



we finally obtain Eq. (|T3D. 



(AlO) 



(All) 



APPENDIX B: EXPLICIT EXPRESSIONS FOR Eq AND Ej 

The energy Eq, calculated from Eq. (p^), contains the divergent term 



47m „ 



G=o 



27ms 27ins ^ ( G 

-G 

27r 



G 



G \2^jfrm's 
— 2e^^fru, 



G=0 



G=0 



where we have made use of the limit 

lim erf(x) = 2/1/7?- 



a:->0 



(Bl) 



(B2) 



in the second line of (|B1|). The divergent term in the last line of ( |B1D is independent of the 
lattice geometry and can be neglected. In fact the divergency is exactly balanced by the 
interaction energy of the electrons with a positive background located in the same layeiH 



dr 



- e n. 



27r 



(B3) 



q=0 



Equation ( ^2]) can also be used to evaluate the contribution to Eq of the last two terms 
in Eq. il) 



lim 

r-+0 



-2y^. 



(B4) 



Using and the identity G = 2TTns (i x R), where z is a unit vector normal to the layers, 
Eq reduces finally to Eq. (]T^), which is equal to Eq. (2.15) of Ref. |1^. Eq/2 gives the static 
energy per electron of a SLWC of density Ug. 
The divergent term in Ej is 



e 



^2 



47rno 



,TTr] 



e 7™, 



G=0 



G 



{ 



eGr?/v^ erf ( — ^ + + e^^^'/v^ erf ( — ^ 

\20m7 j \2y/nns 



(B5) 



G=0 
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The second term on the rhs takes the hmit 



— 2e 7rn, 



TC^/n 



+ 



; erf (v^r/^ 



(B6) 



for G ^ 0. The first term on the rhs of (p5| ) can be rewritten 



G 



+ 



G 



G=0 



G=0 



+ 26^71^/^7/. (B7) 



Again, the divergent term on the rhs is independent of the lattice geometry and can be 
neglected. In fact, this term is exactly balanced by the interaction energy of an electron 
with a positive background charge located at the opposite layer 



dr 



.1/2 



which balances the divergency. Therefore, we obtain 



k 



k=0 



(B8) 



47m o 



G=o 



2e^-yrI7|7rr7erfc (^^/ttt]^ — e '^''^j , 



(B9) 



and, finally, Eq. (P^). 

The background charge does not need to sit in the same layer of the mobile electrons. 
This would be, in fact, the situation in the 2D electron gas realized in semiconductor het- 
erostructures, where the positive ions sit far from the inversion layer. In this case the 
electrostatic energy has an additional contribution 2T[e'^ns{d' + rf"), where d\ d" are the 
distances between the compensating layers and the electron layers. Since this additional 
contribution does not depend on the inter-layer distance d nor on the lattice structure, it 
can be neglected. 
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FIGURES 



FIG. 1. Static energy per particle of phases I, III, and V. In the insets we show the corresponding 
lattice geometries in which dots and crosses identify the two sub-lattices, ai and a2 are given in 
Table |. 

FIG. 2. Detail of Fig. |^ showing the transitions: (a) I^II^III, and (b) III^IV^V. Also shown 
are the lattice geometries in which full dots and crosses identify the two sub-lattices. Empty dots 
and diamonds which indicate the sub-lattices of phases I (a) and III (b) are also reported for 
reference. The insets show how: (a) the aspect ratio a^jax^ and (b) the sine of the angle Q between 
ai and a2 evolve during the transition. 

FIG. 3. Phonon dispersion curves for phase I (top panel), phase III (middle panel), and phase 
V (bottom panel), and for several values of r/, as indicated in the legends. Phonon frequencies are 
shown along high symmetry directions in the Brillouin zones of the three phases. In each panel, 
high symmetry points along the abscissa are labeled according to the insets. Frequencies are given 
in terms of the characteristic frequency = e^n^^"^ /rrie. 

FIG. 4. Sound velocity of the TA mode (o;^ = e^n^/^/mg) along the (1,0) (solid lines) and (1,1) 
(dashed lines) directions for the five phases. The sound velocity of phase V is isotropic and the 
two curves coincide. Vertical dotted lines indicate the phase boundaries, according to Figs. [l| and 
^; the labels on top of the figure indicate which phase is energetically favoured in each region. 

FIG. 5. T=0 phase boundaries (solid dots) and range of stability (crosses) of the five phases 
along the ?7 axis. 

FIG. 6. Optical frequencies at the F point for the five phases. For each phase, LOopt are reported 
in the whole range in which the phase is stable. Vertical dotted lines indicate the phase boundaries; 
the labels on top of the figure indicate which phase is energetically favoured in each region. 

FIG. 7. Phonon DOS as a function of frequency for different values of r]. For each value of rj, 
the DOS corresponding to the energetically favoured lattice is reported. Dashed lines indicate that 
a "soft" phase (either II or IV) is stable at that value of rj. uJl = e^n'^/^/me. 

FIG. 8. Melting temperature Tm for the five phases. For the "soft" phases II and IV (dashed 
lines) we used the value of the continuously changing parameter, either the aspect ratio a^jax or 
the angle between ai and a2, respectively, for which the energy is at its minimum at T=0. 

FIG. 9. Same as Fig. ^, but with Tm plotted in units of {kBd/e^y/T:)~^ . 
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TABLES 



TABLE L Lattice parameters of the five geometries considered, a is the NN distance. For 
each phase, the primitive vectors ai and a2, the inter-lattice displacement c, the reciprocal lattice 
vectors bi and b2, and the charge density Ug are indicated. For phase II, 02/01 is the aspect ratio. 
For phase IV, 9 is the angle between ai and a2. 

Phase ai/o ^2/0 c bi/(27r/a) b2/(27r/a) n^a^ 

I- One-component hexagonal (1,0) (0,^/3) (ai +a2)/2 (1,0) (0,1/^/3) 1/^3 

II- Staggered rectangular (1,0) (0,02/01) (ai-|-a2)/2 (1)0) (0,01/02) 01/02 

III- Staggered square (1,0) (0,1) (ai+a2)/2 (1,0) (0,1) 1 

IV- Staggered rhombic (1)0) (cos ^, sin ^) (ai-|-a2)/2 (1, — cos 0/sin0) (0, l/sin6') 1/sin^ 

V- Staggered hexagonal (1, 0) (1/2, ^/3/2) (ai + aai/S (1,-1/^/3) (0,2/^3) 2/^/3 



TABLE II. Fitting parameters for the sound velocity vta in Eqs.(|3^) and (^). 
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TABLE III. Fitting parameters for the optical frequencies oJopt in Eqs.(^l|), (42), and (|4^). 
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